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The  temperature  distribution  T(x,y,z,t)  in  a  specimen  during 
heating  and  cooling  of  its  outer  surface  is  determined  by  solving  the 
heat  diffusion  equation  for  given  boundary  conditions.  Solutions  are 
often  obtained  for  special  forms,  such  as  infinite  cylinders  of  infi¬ 
nite  slabs,  of  the  specimen,  and  for  some  simplified  boundary  condi¬ 
tions,  such  as  abrupt  initial  temperature  change  at  the  surface  of 
the  sample.  In  practical  problems,  the  solutions  are  often  approxi¬ 
mated  by  simple  forms  of  the  time  temperature  relations.  For 
instance,  in  the  case  of  retorting  of  food,  the  center  temperature 
in  the  can  is  usually  approximated  by  a  zero  order  Bessel  function 
valid  for  very  large  values  of  time.  Such  approximations,  while 
valuable,  are  inadequate  for  exact  studies.  Using  modern  computers, 
we  were  able  to  calculate  T(x,y,z,t)  very  accurately  for  any  practi¬ 
cal  size  cylinder  and  for  different  boundary  conditions  corresponding 
to  Nusselt  numbers  between  0  and  5000.  The  high  accuracy  and  the 
rapid  calculations  make  the  method  very  useful  in  many  fields  of 
thermal  engineering.  In  the  present  paper,  wo  apply  this  method-^-\£  \*\  5  £  4  \  n 
exact  integral  calculations  of  the  survival  fraction  of  bacteria 
during  heat  sterilization  process.  ^ 

2.  HEAT  DIFFUSION  EQUATION  FOR  A  FINITE  CYLINDER 

The  conduction  of  heat  in  a  cylinder  is  given  by  a  second 
order  differential  equation.  In  cylindrical  coordinates  it  is: 

a--^— »  «  k  V2T(r ,  0,  z,  t) 
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where  T(r,  0,  z*  t)  is  the  temperature  at  point  P(r,  0,  z)  and  at 
time  t,  and  ic  is  the  thermal  diffusivity  of  the  material  of  the  solid 

k  -  *  *  thermal  conductivity  (2) 

P*c  density  •  specific  heat 

We  will  consider  heating  a  cylindrical  specimen  from  an 
initial  temperature  with  a  constant  temperature  T^  on  the  outside 
of  the  surface.  Then  the  temperature  T(r*  0,  z*  t)  at  any  point  in 
the  specimen  at  any  time  t  will  be  given  by  Eq.  (1),  together  with 
the  boundary  conditions.  To  simplify  the  equations*  we  introduce  a 
new  temperature  variable: 

T  -  T  -  Th  (3) 

Eq.  (1)  becomes 

f£-  «  ?2T  for  r  <  a,  |z|  <  l  (4) 

where  a  is  the  radius*  and  i  the  half-length  of  the  cylinder.  The 
initial  value  of  T  at  any  point  of  the  specimen  is: 

Ti  -  Ti  -  Th,  for  t  <  o  (5) 

We  will  consider  here  the  case  of  convective  heating  outside 
the  cylinder.  At  the  boundary  there  will  be  an  interface  or  transi¬ 
tion  layer  in  which  the  temperature  is  changing  from  to  the  actual 
temperature  of  the  specimen  at  the  surface.  We  then  obtain  the  fol¬ 
lowing  boundary  conditions  by  equating  the  heat  flowing  through  the 
transition  layer  and  the  heat  flowing  through  the  outermost  layer  of 
the  specimen. 


fl  +  ^T-oatr-a  (6) 

+  (7) 

-|I  +  Jl  T  -  o  at  z  -  (8) 

oZ  A 


where  X  is  conductivity  of  the  specimen*  and  h  *  X gas/ <5 gas,  the  con¬ 
ductivity  of  the  gas  (air)  divided  by  the  thickness  of  the  gas  (air) 
film*  is  the  coefficient  of  heat  transfer  due  to  the  gas  film  only. 
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Due  Co  radiation  at  the  surface,  h  actually  would  be  equal 
to  Agas/6gas  plus  some  constant.  Also,  in  the  case  of  a  meat  sample 
with  casing,  we  have  two  transition  layers  Instead  of  one  at  the  sur¬ 
face,  the  solid  casing  material  and  the  gas  film  adhered  onto  it. 

All  these  three  factors  can  be  combined  into  a  single  coefficient  in 
front  of  T  in  Eqs.  (6),  (7)  and  (8).  Henceforth,  h/A  in  Eqs.  (6), 

(7)  and  (8)  will  stand  for  this  effective  coefficient. 

Eq.  (4)  with  boundary  conditions  Eqs.  (6)  -  (8),  i.e.,  3T/3r 
being  proportional  to  T,  can  be  solved  with  the  usual  technique  of 
the  separation  of  variables. ^ *^)  The  solution  for  the  temperature 
distribution  in  a  finite  cylinder,  in  our  case,  takes  the  following 
form: 

T  »  1  l  Aj,n  *  cos  <XJ*>  *  Jo(«nr>*«““<X5  +  (9) 

J  n 

where  Aj>n  are  the  coefficients  of  expansion  and  where  xn  •  a  •  on 
are  the  roots  of  the  following  equation  (J0  and  Ji  being  Bessel  func¬ 
tions  of  order  zero  and  one  respectively), 

*n  ‘  J1  <*n>  "  f  ‘  Jo  (*n>  (10) 

and  yj  -  Ajt  are  the  roots  of 

M  cos  yj  -  yj  sin  y^  -  o  (11) 


We  have  deduced  from  the  boundary  conditions  the  coefficients 
of  expansion  Aj,n  in  Eq.  (9).  After  some  lengthy  calculation,  Aj  n 
are  found  to  be 


A 


j  ,n 


4  8in  yj  (i2) 


sin(2  yj)  +  2  Yj 


The  temperature  T  -  T(r,  0,  z,  t)  of  the  specimen  is  then 


T  -Th 


2(Ti-Th)  l 
1-1 


J 

o 


<v> 


4  sin  yj 

sin(2  yj)  +  2  yj 


cos  AjZ 


J1  (xn) 

(13) 

e-’f(xj2  ♦  °n2)t 
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where 


v  .  IL  .  a  »  effective  conduction  Nusselt  Number 

*  or  Biot  NuraberO)  (14) 

and  xn  and  are  roots  of  Eqs.  (10)  and  (11). 

3.  CALCULATIONS  OF  THE  TEMPERATURE  DISTRIBUTION 

The  temperature  distribution  T(r,  0,  z,  t)  can  be  computed 
for  any  point  P(r,  0,  z)  inside  the  cylinder  and  any  time  t  provided 
the  roots  xn  of  Eq.  (1)  for  a  given  v,  and  the  roots  y^  of  Eq.  (11) 
for  a  given  value  of  hi/X  are  known. 

Our  task  then  is  to  devise  first  some  means  to  compute 
rapidly  these  roots  with  any  practical  values  of  v  and  hi/X. 

We  have  written  computer  programs  which  enable  us  to  compute 
the  first  36  roots  of  Eqs.  (10)  and  (11),  for  v  =  0  to  5000  and 
hi /X  -  0  to  10,000,  with  rapid  convergence. 

Likewise,  the  double  summation  of  Eq.  (13)  has  been  carried 
out  with  the  computer.  This  double  summation  in  principle  could  be 
carried  out  to  as  many  terms  as  we  wish,  so  long  as  the  values  of  the 
roots  xn  and  y^  are  available. 

We  have  also  used  a  temperature  variable  Y,  the  ,f relative 
excess  temperature** ,  which  is  independent  of  the  actual  initial  and 
heating  temperature , 

*  “  *  I  T.h.  (15) 

V  Th 

It  is  seen  from  Eq.  (13)  that  Y  is  equal  to  twice  the  double 
summation  of  that  equation,  and  is  Independent  of  and  T^.  Thus, 
the  actual  temperature  distributions  T(r,  0,  z,  t) ,  due  to  different 
initial  and  heating  temperatures,  can  be  directly  compared  or  derived 
from  one  another. 

The  temperature  distribution  T(r,  0,  z,  t)  in  specimens  dur¬ 
ing  heating  and  cooling  has  often  been  obtained  by  some  approximations 
for  the  heat  diffusion  equation  or  its  solution.  In  the  present  paper 
we  deduced,  on  the  other  hand,  the  exact  analytical  solution  to  the 
heat  diffusion  equation  for  finite  cylinders  in  convective  heating. 

The  truncation  of  the  two  series  in  Eq.  (13)  to  36  terms  in  each  of 
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the  summations  results  in  a  very  accurate  form,  allowing  T(r,  6,  z,  t) 
to  be  calculated  for  any  practical  size  finite  cylinders  and  for 
Nusselt  (Biot)  number  extending  from  0  to  5000*  The  method  developed 
should  have  great  applicability  in  thermal  engineering  Industry. 

Tables  I  -  II  give  the  relative  excess  temperature  Y  on  the  cross 
sections  at  z/i  -  0,  and  0.6  respectively,  at  r/a  *  0,  0.2,  0.4, — 

1.0  and  times  t  *  5,  10,  15, - 240  minutes,  for  a  specimen  of  radius 

a  -  5  cm,  half-length  l  “  2.5  cm,  with  k  *  1.35  x  10~3  cm2s-1  and 
v  -  6.0.  Similar  tables  for  Y  on  the  cross  sections  z/fc  ■  0.2,  0.4, 
0.8,  and  1.0,  are  available  from  the  authors. 

The  last  columns  of  these  tables  give  the  survival  fraction 
calculations  of  bacteria  during  sterilization  to  be  discussed  in 
later  sections. 

Figs,  la-lf  are  the  computer  plots  of  these  tables. 

4.  CALCULATION  OF  BACTERIA  SURVIVAL  IN  HEAT  STERILIZATION 

We  now  apply  the  above  calculations  of  heat  diffusion  to 
thermal  killing  of  bacteria  in  finite  cylindrical  specimens.  The 
microbial  kill  is  described  quantitatively  by  the  survival  fraction 
N/NQ  where  N  is  the  number  of  bacteria  survived  at  time  t  and  N  is 
the  initial  value  of  N.  For  a  large  class  of  microorganisms,  N/N0 
during  thermal  sterilization  at  a  given  constant  temperature  T  follows 
closely  the  following  equation  (1*5,6). 

in  1!_  -  -C.  [exp  (-Ea/RT) ]  •  t  (16) 

No 

where 

Ea  «  the  Arrhenius  activation  energy  in  kcal-mol~1  for 

inactivating  a  unit  or  a  molecule  that  is  essential 
for  the  survival  of  a  bacterium 

R  *  the  universal  gas  constant  in  J-mol*1 

T  ■  sterilization  temperature  in  Kelvin 

t  *  sterilization  time  in  s 

C  ■  a  constant 
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In  this  case  where  T  is  constant,  we  obtain  by  differentia¬ 
tion  of  Eq.  (16),  that  the  kill  rate  function  (dN/dt)/N  is  the  expon¬ 
ential  exp  (-Ea/RT) ,  that  is: 

M  .  I  =*  -C* [exp(-E  /RT) ]  -  -C*F(T(t)  ,  t)  (17) 

dt  N  a 

In  some  experiments  it  is  found  that  the  semi-logarithmic 
plot  of  the  survival  fraction  has  a  "shoulder".  This  may  be  repre¬ 
sented  by  the  function 

F(T(t),t)  -  1  -  [l-exp(-Ea/RT(t))]n  (18) 

which  for  n  -  1  reduces  to  the  form  in  Eq.  (17),  the  frequently  used 
Arrhenius  function: 

F(T(t) ,  t)  =*  exp  [-Ea/RT(t)]  (19) 

We  now  show  that  for  the  general  case  where  the  kill  rate 
function  assumes  a  general  form  F(T(t)  ,  t)  ,  with  the  temperature  T(t) 
varying  with  time  t,  the  survival  fraction  is  given  by  the  following 
time  integral: 

t 

In  H  -  -C  /  F(T(t),  t)  dt  (20) 

N  ° 

From  the  emperical  law  Eq.  (16)  for  constant  T, 

In  ■  -C  •  [exp(-Ea/RT)]  •  T  (16) 

No 

we  obtain  the  following  differential  form  in  terms  of  the  kill  rate 

-  -C  •  [exp(-Ea/RT) ]  •  N  (21) 

dt 

Now  we  generalize  Eq.  (21)  to  include  the  case  where  T  varies 
with  t.  By  integrating  Eq.  (21)  for  this  temperature  varying  case, 
we  have  the  following  time  integral  for  the  logarithmic  survival 
fraction 

t 

Jin  —  ■  -C  /  exp  (-E./RT)  dt  (22) 

N„  o'  a 
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Next,  we  generalize  Eqs.  (21)  for  the  general  kill  rate 
function  F(T(t) ,  t) , 

M  -  _c  .  F(T(t) ,  t)  •  N 
at 

the  logarithmic  survival  fraction  is  then 

log  5.  -  -c  /  F(T(t)  ,  t)  dt 
No  ° 

which  is  Eq.  (20). 


(21a) 


(22a) 


5.  CALCULATIONS  OF  BACTERIAL  SURVIVAL  EQUATION 

We  now  put  our  function  T(t)  of  Eq.  (13)  into  the  "bacterial 
survival  integral"  of  Eqs.  (22)  or  (22a)  to  calculate  the  survival 
fraction  N/NQ  of  the  bacterium.  In  this  paper,  we  carried  out  the 
calculations  for  Cl.  botulinum  spores  for  which  we  assume  the  kill 
rate  function  to  be  given  by  the  exponential  function  exp(-Ea/RT) 
rather  than  the  general  function  F(T(t) ,  t)  which  may  represent  any 
empirical  function  or  other  theoretical  function  such  as  Eyring’s. 


The  calculations  for  other  microorganisms  will  be  exactly 

the  same. 

The  integrand  of  Eq.  (22)  after  substituting  T  from  Eq.  (13) 
looks  quite  complicated.  We  have  written  computer  programs  to  carry 
out  the  computation.  But,  before  we  carry  out  the  computation,  the 
two  constants,  which  characterize  the  bacterium  being  considered,  C 
and  Ea  in  Eqs.  (16)  and  (22),  have  to  be  determined. 

6.  DETERMINATION  OF  Ea 


We  show  in  some  details  below  the  relationship  between  the 
so-called  z-value  in  the  terminology  used  by  the  microbiologist,  the 
temperature  T,  and  the  Arrhenius  activation  energy  Ea. 

The  reaction  rate  equation,  Eq.  (21a), 

-  -kN  -  -C-  [exp(-E  /RT)]  •  N  (21) 

dt  a 

means  that  the  temperature  dependence  of  k  is  through  T  in  the  expon¬ 
ential  function  exp(-E  /RT)  only  where  Efl  is  a  constant  independent 
of  T.  Otherwise,  k  will  have  to  be  expressed  by  the  general  function 
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F(T(t)  ,  t).  Thus,  at  temperatures  Tj  and  T2,  we  have 

Ol  "  -C'  CexP(-Ea/RTl)]  ■  N 


O2  “  “C'  [exP("Ea/RT2>l  ’  N 

From  Eqs.  (23)  and  (24), 

Cd»/dt)i  E,  T2  -  T 

—  -  exp  I” - ] 


(dN/dt). 


(23) 

(24) 

(25) 


R  TjTj 


Now  if  we  choose  the  temperature  T2  such  faat  the  rate  of 
Inactivation  is  changed  by  a  factor  of  10  from  that  at  Tt,  that  is: 


(dN/dt) ! 


10 


(dN/dt)2 

whence,  using  Eq.  (25)  it  follows  that: 

E  ,  (in  10) *R  TjT2  ^  (in  10) *R  TjT2 

a _ _ 

Ti  -  T2  z 


(26) 


(27) 


where,  by  definition,  -  T2 ,  the  temperature  change  as  specified, 
is  the  z-value. 

In  case  of  Cl.  botulinum  spores,  representative  values  for 
z  are  in  the  range  z  *  10.4  +  0.8°  at  121.1°c(^).  From  Eq.  (27),  we 
have  then  for  Cl.  botulinum  spores. 


E  -  (in  10)  •  394.26  •  384.26  •  1.987  •  10~3 


10.4 

*  66.6+5.4  kcal.  mol-* 
7.  DETERMINATION  OF  THE  CONSTANT  C 


(28) 


For  Cl.  botulinum  spores,  we  will  use  for  12D  (i.e.,  reduc¬ 
tion  of  the  spore  number  to  10“12  of  the  Initial  number)  the  conser¬ 
vative  FQ-value  of  3.5  min  for  non-acid  and  non-cured  meats,  which 
means  that  heating  at  121. 1°C  for  17.5  sec.  (■  3.5  x  60/12)  reduces 
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the  number  N  by  one  order  of  magnitude.  Thus,  from  Eq.  (16),  for 
Cl.  botulinum  spores. 


[C*exp(-E _/RT) ]  •  17.5  -  in  10 
s 


giving 

C  -  2.3026  e  85.076 
17.5 

=  1.17  •  1036s~1 

Calculations  are  made  rather  simply  with  the  computer  of 
Eq.  (22)  as  a  function  of  time.  The  last  columns  of  Tables  I  and  II 
•give  the  values  of  the  survival  fraction  N/N0  at  the  axis  thus  com¬ 
puted.  Fig.  2  is  the  plot  of  these  "integral"  survival  curves  for 
Cl.  botulinum  spores  at  the  center  of  the  cross-sectional  planes 
z fl  *  0,  0.4,  0.6,  0.8  and  1.0,  of  the  cylinder  with  a  *  5  cm,  v  *  6 
and  k  *  1.35  x  10~3  cm2s_1.  The  heating  and  the  initial  temperatures 
T^  and  T^  are  121. 1°C  and  21°C  respectively. 

It  is  seen  from  the  survival  curve  for  z  ■  o  that  in  order 
to  reduce  the  survival  fraction  of  Cl_.  botulinum  spores  to  10~12  of 
the  initial  spore  number,  about  1-1/2  hours  heating  is  needed  for  a 
beef  roll  of  radius  3  5  cm  and  5  cm  long. 
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